O objetivo desse notebook é aplicar a Linguagem R em Ciência de Dados na resolução de problemas de forma não supervisionada.

kMeans

O algoritmo de agrupamento kMeans é um método de quantização vetorial, originário do processamento de sinais, que visa particionar n observações em k clusters em que cada observação pertence ao cluster com a média mais próxima (centros de cluster ou centróide de cluster), servindo como um protótipo de cluster. Isso resulta em um particionamento do espaço de dados.

Execução do kMean (Adaptado da Wikepedia)

O algoritmo kMeans clássico pode ser definido como a soma dos mínimos locais gerados pelos agrupamentos definido como:

Formulação matemática do kMeans

Base USArrests

A base de dados USArrests contém estatísticas por 100.000 habitantes de prisões por agressão, assassinato e estupro em cada um dos 50 estados dos Estados Unidos em 1973.

data('USArrests')
usa = USArrests
head(usa)

É muito importante tratar valores ausêntes. A função is.na nos diz se um elemento é NA. Combinada com a função any, podemos descobrir se algum elemento da base é NA.

head(is.na(usa))
##            Murder Assault UrbanPop  Rape
## Alabama     FALSE   FALSE    FALSE FALSE
## Alaska      FALSE   FALSE    FALSE FALSE
## Arizona     FALSE   FALSE    FALSE FALSE
## Arkansas    FALSE   FALSE    FALSE FALSE
## California  FALSE   FALSE    FALSE FALSE
## Colorado    FALSE   FALSE    FALSE FALSE
print(any(is.na(usa)))
## [1] FALSE

Também é muito importante normalizar os dados. Poderiamos usar normalização por maxmin.

maxmin <- function(x) {
  return ((x - min(x))/(max(x) - min(x)))
}

normalize <- function(df) {
  return (sapply(df, maxmin))
}

usa = normalize(usa)
head(usa)
##         Murder   Assault  UrbanPop      Rape
## [1,] 0.7469880 0.6541096 0.4406780 0.3591731
## [2,] 0.5542169 0.7465753 0.2711864 0.9612403
## [3,] 0.4397590 0.8527397 0.8135593 0.6124031
## [4,] 0.4819277 0.4965753 0.3050847 0.3152455
## [5,] 0.4939759 0.7910959 1.0000000 0.8604651
## [6,] 0.4277108 0.5445205 0.7796610 0.8113695

Mas vamos usar a normalização pela média e desvio padrão.

usa = scale(usa)
head(usa)
##          Murder   Assault   UrbanPop         Rape
## [1,] 1.24256408 0.7828393 -0.5209066 -0.003416473
## [2,] 0.50786248 1.1068225 -1.2117642  2.484202941
## [3,] 0.07163341 1.4788032  0.9989801  1.042878388
## [4,] 0.23234938 0.2308680 -1.0735927 -0.184916602
## [5,] 0.27826823 1.2628144  1.7589234  2.067820292
## [6,] 0.02571456 0.3988593  0.8608085  1.864967207
summary(usa)
##      Murder           Assault           UrbanPop             Rape        
##  Min.   :-1.6044   Min.   :-1.5090   Min.   :-2.31714   Min.   :-1.4874  
##  1st Qu.:-0.8525   1st Qu.:-0.7411   1st Qu.:-0.76271   1st Qu.:-0.6574  
##  Median :-0.1235   Median :-0.1411   Median : 0.03178   Median :-0.1209  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.7949   3rd Qu.: 0.9388   3rd Qu.: 0.84354   3rd Qu.: 0.5277  
##  Max.   : 2.2069   Max.   : 1.9948   Max.   : 1.75892   Max.   : 2.6444

Agora nossos dados estão prontos. Vamos aplicar a função de agrupamento kMeans para 2 grupos.

cls = kmeans(usa, centers=2)
print(cls)
## K-means clustering with 2 clusters of sizes 30, 20
## 
## Cluster means:
##      Murder    Assault   UrbanPop       Rape
## 1 -0.669956 -0.6758849 -0.1317235 -0.5646433
## 2  1.004934  1.0138274  0.1975853  0.8469650
## 
## Clustering vector:
##  [1] 2 2 2 1 2 2 1 1 2 2 1 1 2 1 1 1 1 2 1 2 1 2 1 2 2 1 1 2 1 1 2 2 2 1 1 1 1 1
## [39] 1 2 1 2 2 1 1 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 56.11445 46.74796
##  (between_SS / total_SS =  47.5 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Analisando os grupos:

cls$cluster
##  [1] 2 2 2 1 2 2 1 1 2 2 1 1 2 1 1 1 1 2 1 2 1 2 1 2 2 1 1 2 1 1 2 2 2 1 1 1 1 1
## [39] 1 2 1 2 2 1 1 1 1 1 1 1

Conseguimso descobrir qual o número de grupos ideal? Uma alternativa é utilizar uma medida de validação de agrupamento como a silhueta. Para isso precisamos carregar a biblioteca cluster.

require('cluster')
## Loading required package: cluster
sil = silhouette(cls$cluster, dist(usa))
print(sil)
##       cluster neighbor sil_width
##  [1,]       2        1 0.3354254
##  [2,]       2        1 0.3254396
##  [3,]       2        1 0.3911248
##  [4,]       1        2 0.1345762
##  [5,]       2        1 0.3774168
##  [6,]       2        1 0.2972633
##  [7,]       1        2 0.5306315
##  [8,]       1        2 0.1743774
##  [9,]       2        1 0.5035739
## [10,]       2        1 0.4044931
## [11,]       1        2 0.4194608
## [12,]       1        2 0.5680837
## [13,]       2        1 0.3354701
## [14,]       1        2 0.4428995
## [15,]       1        2 0.5945178
## [16,]       1        2 0.5362530
## [17,]       1        2 0.3496924
## [18,]       2        1 0.4387476
## [19,]       1        2 0.5625723
## [20,]       2        1 0.4890198
## [21,]       1        2 0.3837284
## [22,]       2        1 0.5165960
## [23,]       1        2 0.5986325
## [24,]       2        1 0.2620747
## [25,]       2        1 0.1134541
## [26,]       1        2 0.5240161
## [27,]       1        2 0.5926288
## [28,]       2        1 0.4299099
## [29,]       1        2 0.5870108
## [30,]       1        2 0.1782450
## [31,]       2        1 0.5222275
## [32,]       2        1 0.3865530
## [33,]       2        1 0.2596933
## [34,]       1        2 0.5143664
## [35,]       1        2 0.3618181
## [36,]       1        2 0.4064074
## [37,]       1        2 0.1942476
## [38,]       1        2 0.5253841
## [39,]       1        2 0.3635329
## [40,]       2        1 0.3650274
## [41,]       1        2 0.5325247
## [42,]       2        1 0.3175580
## [43,]       2        1 0.3471312
## [44,]       1        2 0.4078063
## [45,]       1        2 0.4408083
## [46,]       1        2 0.2622062
## [47,]       1        2 0.3313743
## [48,]       1        2 0.4592799
## [49,]       1        2 0.5891363
## [50,]       1        2 0.4400334
## attr(,"Ordered")
## [1] FALSE
## attr(,"call")
## silhouette.default(x = cls$cluster, dist = dist(usa))
## attr(,"class")
## [1] "silhouette"

Podemos mostrar graficamente o valor da silhueta por amostra.

plot(sil)

Ou usar a média para avaliar se o valor é bom.

mean(sil[,3])
## [1] 0.408489

Podemos variar o número de grupos para tentar encontrar o valor ótimo de grupos baseado na medida silhueta.

cls = sapply(2:10, function(k) {
  aux = kmeans(usa, centers=k)
  sil = silhouette(aux$cluster, dist(usa))
  mean(sil[,3])
})

Podemos plotar os valores da silhueta.

names(cls) = 2:10
plot(names(cls), cls, type = "l")

Usando a biblioteca plotly.

require(plotly)
## Loading required package: plotly
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
data = data.frame(cls, k=2:10)
plot_ly(data, x=~k, y=~cls, type='scatter', mode='lines')

Vamos carregar duas bibliotecas para nos ajudar a visualizar os dados.

require('gridExtra')
## Loading required package: gridExtra
require('factoextra')
## Loading required package: factoextra
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Vamos visualizar os agrupamentos.

plot = lapply(2:5, function(k) {
  cls = kmeans(usa, centers=k)
  fviz_cluster(cls, data=usa)
})
do.call("grid.arrange", plot)

Comprimindo Áudio com kMeans

Um experimento interessante é a tarefa de comprimir áudio utilizando o algoritmo de agrupamento kMeans. Para isso vamos precisar das bibliotecas audio e seewave.

require('audio')
## Loading required package: audio
require('seewave')
## Loading required package: seewave
## 
## Attaching package: 'seewave'
## The following object is masked from 'package:plotly':
## 
##     export

Vamos baixar uma música da internet com permissão de reprodução pelo pixabay. Essa música tem o formato MP3.

url = 'https://cdn.pixabay.com/download/audio/2021/04/07/audio_0bed53371b.mp3?filename=blues-vibes-100-bpm-michael-kobrin-3780.mp3'
download.file(url, 'audio.mp3')

Na sequencia, precisamos converter para WAV e pegar uma amostra para podermos demostrar a compressão. Para isso vamos usar o programa ffmpeg.

system('ffmpeg -i audio.mp3 audio.wav')
system('ffmpeg -i audio.wav -ac 1 mono.wav')
system('ffmpeg -i mono.wav -ss 00:00 -to 00:10 out.wav')

Podemos carregar a música.

sound = load.wave("out.wav")
class(sound)
## [1] "audioSample"

E plotar…

plot(sound, type = "l")

Usando plotly…

data = data.frame(as.numeric(sound), x=1:length(sound))
plot_ly(data, x=~x, y=~sound, type='scatter', mode='lines')

Verificando quantos valores únicos temos nessa amostra.

length(unique(as.numeric(sound)))
## [1] 44923

Para aplicar o kMeans, precisamos transformar esses dados.

wave = data.frame(amp=as.numeric(sound))
summary(wave)
##       amp            
##  Min.   :-0.9010590  
##  1st Qu.:-0.0736106  
##  Median : 0.0000000  
##  Mean   :-0.0000319  
##  3rd Qu.: 0.0752258  
##  Max.   : 0.9559326
head(wave)

Agora podemos aplicar o kMeans. Vamos tentar agrupar os sinais de amplitude da música de 10 até 100 grupos.

require('parallel')
## Loading required package: parallel
cls = mclapply(seq(10, 100, by=10), mc.cores=6, function(k) {
  kmeans(wave, centers=k)
})

O resultado é algo como:

cls[[1]]$centers
##             amp
## 1   0.404090231
## 2  -0.442804785
## 3   0.237218137
## 4  -0.292478713
## 5   0.646721157
## 6  -0.074450045
## 7   0.004286873
## 8   0.112120472
## 9  -0.178236815
## 10 -0.664722312
head(cls[[1]]$cluster, 100)
##   [1] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
##  [38] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
##  [75] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
cmp = cls[[1]]$centers[cls[[1]]$cluster]
plot(cmp, type = "l")

Agora, podemos salvar alguns desses resultados do kMeans para compressão de música para diferentes valores de k e testar.

cmp = cls[[1]]$centers[cls[[1]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk10.wav')
cmp = cls[[2]]$centers[cls[[2]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk20.wav')
cmp = cls[[5]]$centers[cls[[5]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk50.wav')
cmp = cls[[10]]$centers[cls[[10]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk100.wav')

Comprimindo Imagens com kMeans

Um experimento interessante é a tarefa de comprimir imagens utilizando o algoritmo de agrupamento kMeans. Para a compressão, o número de cores precisa ser reduzido. A redução é feita através do agrupamento das cores e seleção daquelas mais frequentes. Para isso precisamos carregar a blbioteca imager:

require('imager')
## Loading required package: imager
## Loading required package: magrittr
## 
## Attaching package: 'imager'
## The following object is masked from 'package:magrittr':
## 
##     add
## The following object is masked from 'package:audio':
## 
##     play
## The following object is masked from 'package:plotly':
## 
##     highlight
## The following objects are masked from 'package:stats':
## 
##     convolve, spectrum
## The following object is masked from 'package:graphics':
## 
##     frame
## The following object is masked from 'package:base':
## 
##     save.image

Algumas imagens são pré-carregadas pelo pacote como a imagem de barcos:

plot(boats, axes=FALSE)

Vamos salvar a imagem original em um arquivo e imprimir seu tamanho em disco:

save(boats, file='boats1.zip')
cat('Tamanho da imagem original:', file.info("boats1.zip")$size, 'bytes')
## Tamanho da imagem original: 1499621 bytes

Podemos extrair informações sobre o objeto da seguinte forma:

boats
## Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 3
dim(boats)
## [1] 256 384   1   3

Gerando o data.frame da imagem boats:

image = as.data.frame(boats)
# cabecalho do data.frame
head(image)

Podemos deixar a imagem borrada:

plot(isoblur(boats, 10), axes=FALSE)

Podemos deixar a imagem em escala de cinza:

plot(grayscale(boats), axes=FALSE)

Podemos remover a cor vermelha da imagem:

test_image = image
test_image[test_image$cc == 1, 'value'] = 0
test_image = as.cimg(test_image)
## Warning in as.cimg.data.frame(test_image): Guessing image dimensions from
## maximum coordinate values
plot(test_image, axes=FALSE)

Podemos separar as cores RGB da imagem da seguinte forma:

red = image[image$cc == 1, 'value']
blue = image[image$cc == 2, 'value']
green = image[image$cc == 3, 'value']

Reconstruindo a base de dados em um data.frame:

data = data.frame(red, blue, green)

Aplicando o kMeans com 10 centroides e 500 iterações máximas:

centroides = kmeans(data, centers=10)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)
summary(centroides)
##              Length Class  Mode   
## cluster      98304  -none- numeric
## centers         30  -none- numeric
## totss            1  -none- numeric
## withinss        10  -none- numeric
## tot.withinss     1  -none- numeric
## betweenss        1  -none- numeric
## size            10  -none- numeric
## iter             1  -none- numeric
## ifault           1  -none- numeric

Adicionando os centroides encontrados em data:

data[,"cluster"] = centroides$cluster
head(data)

O próximo passo é reconstruir a imagem utilizando as cores identificadas como centroides. A imagem será reconstrída com as 10 cores principais:

new_data = do.call('rbind',
  lapply(1:nrow(data), function(i) {
    centroides$centers[data[i,4],]
  })
)

Plotando a nova imagem:

new_image = image
new_image$value = c(new_data[,1], new_data[,2], new_data[,3])
new_image = as.cimg(new_image)
## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values
plot(new_image, axes=FALSE)

Vamos salvar a nova imagem em um arquivo e imprimir o seu tamanho comprimido em disco:

save(new_image, file='boats2.zip')
cat('Tamanho da imagem comprimida:', file.info("boats2.zip")$size, 'bytes')
## Tamanho da imagem comprimida: 61178 bytes

Antes de continuar com o experimento, vamos construir a função imageGenerator que recebe uma imagem e o número de centroides e retorna uma nova imagem comprimida:

imageGenerator <- function(image, centers) {
  
  image = as.data.frame(image)
  red = image[image$cc == 1, 'value']
  blue = image[image$cc == 2, 'value']
  green = image[image$cc == 3, 'value']
  data = data.frame(red, blue, green)
  
  centroides = kmeans(data, centers)
  data[,"cluster"] = centroides$cluster

  new_data = do.call('rbind',
    lapply(1:nrow(data), function(i) {
      centroides$centers[data[i,4],]
    })
  )
  
  new_image = image
  new_image$value = c(new_data[,1], new_data[,2], new_data[,3])
  as.cimg(new_image)
}

Usando a função imageGenerator podemos gerar diferentes imagens com diferentes números de cores. Segue o resultado do experimento para 10, 32, 64 e 128 cores:

par(mfrow=c(2,2))
plot = sapply(c(10, 32, 64, 128), function(i) {
  plot(imageGenerator(boats, i), axes=FALSE)
})
## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)
## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)
## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values
## Warning: did not converge in 10 iterations
## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values